Finite-size scaling in globally coupled phase oscillators with a general coupling scheme 
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We investigate a critical exponent related to synchronization transition in globally coupled non- 
identical phase oscillators. The critical exponents of susceptibility, correlation time, and correlation 
size are significant quantities to characterize fluctuations in coupled oscillator systems of large but 
finite size and understand a universal property of synchronization. These exponents have been 
identified for the sinusoidal coupling but not fully studied for other coupling schemes. Herein, for a 
general coupling function including a negative second harmonic term in addition to the sinusoidal 
term, we numerically estimate the critical exponent of the correlation size, denoted by pj^, in a 
synchronized regime of the system by employing a non-conventional statistical quantity. First, we 
confirm that the estimated value of v+ is approximately 5/2 for the sinusoidal coupling case, which 
is consistent with the well-known theoretical result. Second, we show that the value of i/+ increases 
with an increase in the strength of the second harmonic term. Our result implies that the critical 
exponent characterizing synchronization transition largely depends on the coupling function. 

PACS numbers: 64.60.F-, 05.45.Xt 
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I. INTRODUCTION 

Populations of coupled rhythmic elements can exhibit synchronization and collective behavior via mutual interac- 
tions [ij. Such phenomena are observed in a variety of systems such as chemical reactions, engineering circuits, and 
biological populations jM- To elucidate the general properties of such phenomena, the phase description of systems 
has been widely used In particular, there have been many studies on globally coupled phase oscillators 0,01, 

defined as follows: 
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where 9j is the phase of the jth oscillator, is the natural frequency of the jth oscillator, X > is the coupling 
strength, h is the coupling function, and N is the number of oscillators. When h{x) — sin(a:), this model is referred to 
as the Kuramoto model j^] . One of the main issues in this model has been the scaling property of the order parameter 
defined as follows 0: 
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In the thermodynamic limit {N — > oo), the phase oscillator model ([T|) exhibits a synchronization transition when 
the coupling strength K surpasses a critical value Kc- This transition can be characterized by a change of the order 
parameter from zero to a non-zero value. We assume that the stationary state (i?(t) — 0) in the incoherent regime 
supercritically bifurcates at the critical coupling strength K = Kc, above which the oscillators are synchronized. 
The behavior of the order parameter is exemplified for a finite-size system in Fig. [T] . The scaling law of the 
order parameter with respect to a change in the coupling strength K around the synchronization transition point has 
been well studied in relation to the second order phase transition 0, |3) S 043 ■ The scaling property has been fully 
understood in the case where A'^ is infinite, not only for the sinusoidal coupling function but also for general coupling 
functions 0, However, it is less clear in finite-size systems because the property of fluctuations depends on the 

system size. 
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FIG. 1. Bifurcation diagram of the order parameter 7?, where h(x) 
bifurcation occurs. 



sin(2;) — 0.5 sin(2a;) and A'^ = 8000. A supercritical 



The critical exponents of the order parameter, correlation time, susceptibility (corresponding to the product of the 
variance of the order parameter and A^) , and correlation size (representing the number of oscillators which are almost 
synchronized but not completely) are significant quantities that can characterize the behavior near phase transitions 
in physical systems. The values of the critical exponents of these statistical quantities can be used to categorize 
physical systems into minor classes, because the values are independent of the details of the system [l^]- Further, in 
equilibrium systems, the critical exponents of these statistical quantities completely determine those of all the other 
statistical quantities ■ Therefore, the evaluation of the critical exponents in the phase oscillator model ([T]) enables 
the clarification of differences between equilibrium and non-equilibrium systems. In particular, fluctuations in the 
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systems of large but finite size can be characterized by the exponents of susceptibility, correlation time, and correlation 
size. Although the critical ex pon ents in the phase oscillator model (IT|) with finite large N have been obtained for the 
sinusoidal coupling function [lll - ll6j . it is known that the critical exponent of the order parameter depends on the 
coupling scheme 0, • This fact motivated us to examine if the critical exponents of other statistical quantities also 
depend on the coupling function or not. 

In the present paper, we employ a non-conventional statistical quantity to evaluate the critical exponent of corre- 
lation size, 1/+, in the synchronized regime of the phase oscillator model ^ with finite larg e N. This is because it is 
difficult to compute the value of iy+ using the critical exponent of the order parameter |16| . The statistical quantity 
that we use is denoted by D, which is given by the diffusion coefficient of the temporal integration of the order pa- 
rameter, multiplied by system size N jl7|. Using the statistical quantity D, we perform the finite-size scaling analysis 
[lOi. First, we confirm that the estimated value of is ap pro ximately 5/2 for the sinusoidal coupling h{x) = sin(a;), 
which is consistent with the well-known theoretical result |l3l [Tsj . Second, we consider a general coupling function 
including a negative second harmonic term in addition to the sinusoidal term, i.e. h{x) = sin(a;) — q sin(2a;) with q > 0. 
We show that the value of increases with an increase in the strength q of the second harmonic term. Our result 
means that the critical exponent characterizing synchronization transition largely depends on the coupling function. 

Although fluctuations of the order parameter in the phase oscillator model ([T]) have been studied for the past two 
decades 11 17], those for a general coupling function have not been fully understood. In particular, for any general 
coupling scheme other than the sinusoidal one, the critical exponents of statistical quantities have not been reported 
except for that of the non-conventional statistical quantity D [l3]. Note that the values of critical exponents are 
independent of the details of systems and they characterize universal structures of the systems [l0|. In fact, in the 
thermodynamic limit — )■ oo, it was previously shown that as far as the coupling function includes a second harmonic 
term, the critical exponent of the order parameter takes the same value [aQ. In contrast, our finding means that 
the value of crucially depends on the strength of the second harmonic term. Our result highlights a new relation 
between the couping schemes and universal structures in the globally coupled phase oscillators ([1]). 



II. STATISTICAL QUANTITY FOR THE FINITE-SIZE SCALING ANALYSIS 
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FIG. 2. The time evolutions of the variance cr^(t) of the integrated order parameter in the coherent regime of the Kuramoto 
model where h{x) — sin(a;), A'^ = 1000, and K = 1.635 > Kc ~ 1.59 • ■ • . These two evolutions are generated for different initial 
conditions. The value of a^{t) increases linearly with slope 2Dt after a transient period in both cases. Because of multistability, 
the value of D differs depending on the initial conditions. 



We consider the diffusion coefficient of the time integral of R{t), which characterizes long-term fiuctuations in R{t) 
[l3]. The variance of R{s)ds is given by: 
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a'{t) = N lim- / Ris)ds~{R)tt] dto, (3) 
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FIG. 3. Probability density function of D estimated by using 300 different initial conditions in the model ([T]), where h(x) = sin(x) 
and K = 1.635 > Kc — 1.59 ■ ■ • . The horizontal axis is normalized by {D)e, where (■)e is the ensemble average over 300 different 
initial conditions. We numerically obtained 300 values of D for each A'^ and estimated the probability density function of D by 
means of kernel density estimation. The standard deviation of D is sufBciently small for each A''. 

where {R)t represents the time average, defined as follows: 

{R)t^\im ^ r R{t)dt. (4) 

The following diffusion law then holds: 

D = lim a^{t)/2t. (5) 

The value of D characterizes how /J R{s)ds differs from its mean value {R)tt on the ensemble average. The value of D 
is estimated by fitting cr^(t) with a line of slope 2Dt except for a transient period as shown in Fig. [2] Different initial 
conditions may yield different values of D because of multistability 191], as seen in Fig. [51 However, the dependence 
of D on initial conditions is sufficiently small as shown in Fig. |21 This result implies that the property on fluctuations 
is almost similar for the multiple coexisting solutions. All of them simultaneously transit from a non-synchronous 
state to a synchronous one as the coupling strength K exceeds the critical value K = Kc for a large system size N 

mm- 

III. THE FINITE-SIZE SCALING ANALYSIS 

This study addresses the coupling function in the form h{x) — sin(a::) — (7sin(2a;) with < g < 0.1. We examine 
this coupling function because it is considered to be a "generic" coupling function H, Q ; provided that the coupling 
function includes the second harmonic term in addition to the sinusoidal one, the critical exponent of the order 
parameter does not depend on other higher-order terms [Sj, Q . 

The scaling hypothesis [l3| states that any quantity A, which shows critical divergence at K = Kc, is scaled for 
K > Kc as follows: 

N-'-+/''+A^^{{K - Kc)N^/''+), (6) 

where and r_|_ represent the critical exponents of the correlation size and A in the synchronized regirne, respectively. 
The function 5" is approximated as "^{x) ~ ux~'^+ + x'^'^ with s > r_|_ and a constant u for large x so that the 
orders with respect to system size N in both sides of Eq. © are consistent. It was shown that for a general coupling 
function h{x) in the phase oscillator model ([l}, the critical exponent of D in the synchronized regime is equal to 
[l3]. Therefore, we replace A by D and then substitute r+ = into Eq. ©. Furthermore, we assume u = in the 
above approximation form of ^!{x) because ^(x) — > as a; ^ oo. As a result, we obtain 
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FIG. 4. Estimation of v+ in the Kuramoto model, (a) [K — Kc) vs. (-D)c, where means the ensemble average over 100 



overlapping samples, (b) {K — Kc)N^^'^''^^ vs. {D)c in a log-log plot, which implies i)^^""' 
A'' = 1000 were not used for the line fitting. 
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FIG. 5. Estimation of in the phase oscillator model ((T| with h{x) = sin(a;) — 0.1sin(2a::). (a) {K — Kc) vs. (-D)e, where 
{■)e means the ensemble average over 100 overlapping samples, (b) {K — Kc)N^^^'^^ vs. {D)e in a log-log plot, which implies 
£,0-00) ^ 2.95. The two left endpoints for = 1000 were not used for the line fitting. 



which means that if the finite-size effect is not so strong, (i) D decreases in a power law fashion with system size N for 
a fixed value of K, and (ii) the numerically computed values of D should fall on a straight line against {K — Kc)N-^/'^+ 
in a log- log plot. 

In numerical simulations, the natural frequencies ujj of the individual oscillators are chosen to satisfy 

j/{N+l)= I ' g{Q)dQ, (8) 

where g{uj) is the Gaussian distribution with mean zero and variance one. To avoid a situation in which the finite-size 
effect yields an erroneous value for we limit the range of K to {Kc <)K^ < K < K+. We select the lower value 
K- such that D decreases in a power law fashion with system size N for K > K^. Furthermore, we choose the 
upper value such that the power law decay of D is kept within the range K- < K < Kj^ in order to avoid a 
large variation in the estimated values of v^. In addition, we use a bootstrap method. First, we calculate the value 
of D by using 100 different initial conditions for each pair of {N, K). We randomly choose M overlapping samples of 

D values from the 100 simulation results and average them. Next, we estimate the value of v^^"^ that minimizes the 
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FIG. 6. The standard deviation of the values of estimated from M samples, where 5 < M < 100. (a) The case of g = 0. 
(b) The case of g = 0.1. In both cases, the standard deviation is very small for a sufficiently large M. 
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FIG. 7. Relationship between the estimated exponent z/^*^' and the strength q of the negative second harmonic term in the 
phase oscillator model ([1} with h{x) = sin(a;) — gsin(2a;). (a) The crosses indicate the mean values of C'^^''"''\ The circles and 
the error bars indicate the averages and the standard deviations of the distributions of respectively, (b) Distributions of 
i>9^' for each q. 



mean squared error between the mean values of D and the fitting line in a log-log plot, where v^"^ represents the 
estimated value of 1^4- using the averaged values of D over the M samples. By repeating this procedure 10000 times, 
we obtain the estimated values of the mean and the standard deviation of A 



(M) 



First, let us consider the Kuramoto model, i.e. q = 0. The estimated mean value of v^^^"^ is given as z)^^""^ 



2.49 



by using the method explained above. Figure \M.a) shows the value of D averaged over 100 overlapping samples for 
each pair of {N, K). The same data points shown in Fig. IDJb) are plotted against [K — Kc)N'^^'^''^'^ m a log-log plot. 

The data are fitted well by the straight line. The estimated exponent i>^^'^°-' ~ 2.49 is consistent with the analytical 
results v+ = 5/2 [U, [l3- Therefore, our method can estimate the exact value of v+ well. 

Next, we investigate the model ([1} with a more general coupling function with q = 0.1. Figure [SJ a) shows the 
averaged values of D over 100 overlapping samples. The estimated mean value of i/^:^'^^^ is given as j>^^*'°' ~ 2.95, as 
seen in Fig. [S^b). This means that the value of in this model is larger than that in the Kuramoto model. 

Note that M = 100 samples are enough to estimate the accurate value of because the standard deviation of the 
estimated values of i>^^^^ is sufficiently small, as shown in FigslHlJa) andlHl^b). 

Finally, we examine the relationship between the exponent and the strength q of the negative second harmonic 
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FIG. 8. Relationship between the diffusion coefficient D and the system size A'^. D takes its maximal value in the coherent 
regime due to the finite-size effect. The value of the coupling strength K exhibiting the maximal value of D gradually approaches 
K = Kc as A'^ increases. A fixed initial condition is used to compute the value of D for each A'^. 



term. The mean of is likely to increase with an increase in the value of q as shown in Fig. [7]Ja). In addition, 

we obtain the distributions of as shown in Fig. [Tfb). The standard deviation of is not so large and the mean 

value of i)^^^ is also likely to increase with an increase in the value of q. These results imply that the critical exponent 
depends on the second harmonic term of the coupling function. 



IV. DISCUSSION 



We explain why our method yields a good estimation for the critical exponent i/j^. In the limit N ^ oo, D ^ oo 
as K ^ Kc — 0, whereas D — for K > Kc jl3|- However, for a finite A^, D takes its maximal value in the coherent 
regime as shown in Fig. |S1 It means that D virtually reflects the feature of the incoherent solution of the infinite-size 
system even for K > Kc if K ~ Kc and A^ is finite. Such a finite-size effect is well known in the literature As a 
result, to well estimate the value of by the finite-size scaling analysis, we must avoid the coherent regime that is 
very close to the transition point K — Kc- Our numerical simulations have selected the region in which the finite-size 
effect is weak enough so that D ^ 0{1/N°'). Note that it is difficult to determine whether the finite-size effect is 
sufficiently weak by using other statistical quantities. 

We remark on another numerical study about the critical exponent of correlation size (Tg} . Let us denote the critical 
exponent of the order parameter i? by /? and that of correlation size in the incoherent regime by u- . It is known that 
R ^ N-^/" for K ^ Keif v+ = i^-{=: v). The value of v can be computed as = 5/4 for the Kuramoto model 
with deterministically chosen natural frequencies Iff]. However, there is little evidence for = v-. In fact, as we 
have already mentioned, the value of the critical exponent of D differs depending on whether the system behavior is 
coherent or incoherent [17| . 

Our numerical simulations have shown that when the coupling function possesses a negative second harmonic term 
with strength q in addition to the sinusoidal one, > 5/2 and increases with q. This result is consistent with 
the following property. Near the synchronization transition point K = Kc^ R is scaled as i? ^ (1-1- l/q){K — Kc) + 
0{{K - Kef) [1. This implies that the more the value of q increases, the more slowly the fiuctuations of this system 
decay with the coupling strength K . 

Our method has the following potential application. The method in this paper enables us to obtain an approximate 
value of for a general coupling function. The value of v+ is the same for finite-size scaling analysis of other 
statistical quantities [l^]- Suppose that we analytically obtain the value of v+ or the critical exponents of other 
statistical quantities such as susceptibility and correlation time for a general coupling function; our methods make it 
possible to confirm the validity of the analytical results numerically by the finite-size scaling analysis [lo| because we 
can compute an approximate value of . 
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